Evolutionary rescue experiment: effect of genetic diversity
Laure Olazcuaga
04-February-2022
1 DATA
1.1 Pop sive over time dataset
data <- read.table(file=here::here("data", "BE_Census_Data.csv"), sep=",", header=TRUE)
head(data)## Block Old_Label Label Population Genetic_Diversity
## 1 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High1 High1 High
## 2 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High2 High2 High
## 3 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High3 High3 High
## 4 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High4 High4 High
## 5 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High5 High5 High
## 6 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High6 High6 High
## Generation_Parents Generation_Eggs Date_Census Date_Start_Eggs Date_End_Eggs
## 1 0 1 2/17/21 2/17/21 2/18/21
## 2 0 1 2/17/21 2/17/21 2/18/21
## 3 0 1 2/17/21 2/17/21 2/18/21
## 4 0 1 2/17/21 2/17/21 2/18/21
## 5 0 1 2/17/21 2/17/21 2/18/21
## 6 0 1 2/17/21 2/17/21 2/18/21
## Image_Box1 Image_Box2 MC_Box1 MC_Box2 MC_Total_Adults HC_Box1 HC_Box2
## 1 DSC_0155 <NA> NA NA 101.9 NA NA
## 2 DSC_0156 <NA> NA NA 104.7 NA NA
## 3 DSC_0157 <NA> NA NA 105.8 NA NA
## 4 DSC_0158 <NA> NA NA 101.9 NA NA
## 5 DSC_0159 <NA> NA NA 102.2 NA NA
## 6 DSC_0160 <NA> NA NA 101.7 NA NA
## HC_Total_Adults Nb_parents Growth_rate First_Throw Extinction_finalstatus
## 1 100 NA NA no no
## 2 100 NA NA no no
## 3 100 NA NA no no
## 4 100 NA NA no no
## 5 100 NA NA no no
## 6 100 NA NA no no
## Less_than_5
## 1 no
## 2 no
## 3 no
## 4 no
## 5 no
## 6 notail(data)## Block Old_Label Label Population
## 571 Block5 B5 - Backup Fam H 8/25 BE B5 | G6 Med 19 Med19
## 572 Block5 B5 - Backup Fam I 8/25 BE B5 | G6 Med 20 Med20
## 573 Block5 B5-Backup Fam A 8/25 BE B5 | G6 Med 21 Med21
## 574 Block5 B5-Backup Fam B 8/25 BE B5 | G6 Med 22 Med22
## 575 Block5 B5-Backup Fam C 8/25 BE B5 | G6 Med 23 Med23
## 576 Block5 B5-Backup Fam D 8/25 BE B5 | G6 Med 24 Med24
## Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 571 Med 5 6 8/25/21
## 572 Med 5 6 8/25/21
## 573 Med 5 6 8/25/21
## 574 Med 5 6 8/25/21
## 575 Med 5 6 8/25/21
## 576 Med 5 6 8/25/21
## Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2 MC_Box1 MC_Box2
## 571 8/25/21 8/26/21 DSC_0273 DSC_0274 2.00000 2.00000
## 572 8/25/21 8/26/21 <NA> <NA> NA NA
## 573 8/25/21 8/26/21 DSC_0277 DSC_0278 29.76978 32.92378
## 574 8/25/21 8/26/21 <NA> <NA> NA NA
## 575 8/25/21 8/26/21 DSC_0267 DSC_0268 25.31021 46.10092
## 576 8/25/21 8/26/21 DSC_0283 DSC_0284 71.09684 55.37443
## MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 571 4.0 2 2 4 8 0.500000
## 572 0.0 NA NA NA 0 NA
## 573 62.7 23 23 46 35 1.314286
## 574 0.0 NA NA NA NA NA
## 575 71.4 20 36 56 40 1.400000
## 576 126.5 76 57 133 41 3.243902
## First_Throw Extinction_finalstatus Less_than_5
## 571 no no yes
## 572 extinct yes yes
## 573 no no no
## 574 extinct yes yes
## 575 no no no
## 576 no no nodim(data)## [1] 576 23summary(data)## Block Old_Label Label Population
## Length:576 Length:576 Length:576 Length:576
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## Length:576 Min. :0.0 Min. :1.0 Length:576
## Class :character 1st Qu.:1.0 1st Qu.:2.0 Class :character
## Mode :character Median :2.5 Median :3.5 Mode :character
## Mean :2.5 Mean :3.5
## 3rd Qu.:4.0 3rd Qu.:5.0
## Max. :5.0 Max. :6.0
##
## Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2
## Length:576 Length:576 Length:576 Length:576
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## MC_Box1 MC_Box2 MC_Total_Adults HC_Box1
## Min. : 0.00 Min. : 0.00 Min. : 0.00 Min. : 0.00
## 1st Qu.: 26.52 1st Qu.: 25.99 1st Qu.: 51.05 1st Qu.: 22.50
## Median : 58.87 Median : 54.29 Median :101.50 Median : 57.00
## Mean : 75.88 Mean : 70.57 Mean :128.37 Mean : 72.23
## 3rd Qu.:106.93 3rd Qu.: 99.40 3rd Qu.:170.25 3rd Qu.:101.00
## Max. :327.40 Max. :299.00 Max. :514.40 Max. :288.00
## NA's :142 NA's :141 NA's :5 NA's :145
## HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## Min. : 0.00 Min. : 0.0 Min. : 0.0 Min. :0.0000
## 1st Qu.: 25.00 1st Qu.: 63.0 1st Qu.: 68.0 1st Qu.:0.4379
## Median : 52.00 Median :100.0 Median :100.0 Median :0.9439
## Mean : 68.08 Mean :133.1 Mean :138.6 Mean :1.4639
## 3rd Qu.: 96.75 3rd Qu.:169.0 3rd Qu.:189.0 3rd Qu.:2.3700
## Max. :290.00 Max. :499.0 Max. :499.0 Max. :6.8571
## NA's :142 NA's :48 NA's :123 NA's :147
## First_Throw Extinction_finalstatus Less_than_5
## Length:576 Length:576 Length:576
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## #Remove populations killed due to the pathogen
data_kill<-data[data$Extinction_finalstatus=="kill",]
dim(data_kill)## [1] 72 23data<-data[data$Extinction_finalstatus!="kill",]
#Add pop growth rate
data$Nb_adults <- data$HC_Total_Adults
data$Lambda <- data$Nb_adults / data$Nb_parents
data$obs<-as.factor(1:nrow(data))
#summary(data)
#Remove
#Check variables
str(data)## 'data.frame': 504 obs. of 26 variables:
## $ Block : chr "Block3" "Block3" "Block3" "Block3" ...
## $ Old_Label : chr "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " ...
## $ Label : chr "BE B3 2/17 | G1 High1" "BE B3 2/17 | G1 High2" "BE B3 2/17 | G1 High3" "BE B3 2/17 | G1 High4" ...
## $ Population : chr "High1" "High2" "High3" "High4" ...
## $ Genetic_Diversity : chr "High" "High" "High" "High" ...
## $ Generation_Parents : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Generation_Eggs : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Date_Census : chr "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
## $ Date_Start_Eggs : chr "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
## $ Date_End_Eggs : chr "2/18/21" "2/18/21" "2/18/21" "2/18/21" ...
## $ Image_Box1 : chr "DSC_0155" "DSC_0156" "DSC_0157" "DSC_0158" ...
## $ Image_Box2 : chr NA NA NA NA ...
## $ MC_Box1 : num NA NA NA NA NA NA NA NA NA NA ...
## $ MC_Box2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ MC_Total_Adults : num 102 105 106 102 102 ...
## $ HC_Box1 : int NA NA NA NA NA NA NA NA NA NA ...
## $ HC_Box2 : int NA NA NA NA NA NA NA NA NA NA ...
## $ HC_Total_Adults : int 100 100 100 100 100 100 100 100 100 100 ...
## $ Nb_parents : num NA NA NA NA NA NA NA NA NA NA ...
## $ Growth_rate : num NA NA NA NA NA NA NA NA NA NA ...
## $ First_Throw : chr "no" "no" "no" "no" ...
## $ Extinction_finalstatus: chr "no" "no" "no" "no" ...
## $ Less_than_5 : chr "no" "no" "no" "no" ...
## $ Nb_adults : int 100 100 100 100 100 100 100 100 100 100 ...
## $ Lambda : num NA NA NA NA NA NA NA NA NA NA ...
## $ obs : Factor w/ 504 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...data$Genetic_Diversity <- as.factor(data$Genetic_Diversity)
data$Block <- as.factor(data$Block)
data$Population <- as.factor(data$Population)
data$Extinction_finalstatus <- as.factor(data$Extinction_finalstatus)
#Order levels
data$Genetic_Diversity <- plyr::revalue(data$Genetic_Diversity, c("Med"="Medium"))
data$Genetic_Diversity <- factor(data$Genetic_Diversity, levels = c("High", "Medium", "Low"))
levels(data$Genetic_Diversity)## [1] "High" "Medium" "Low"data<-droplevels(data) #remove previous levels
str(data)## 'data.frame': 504 obs. of 26 variables:
## $ Block : Factor w/ 3 levels "Block3","Block4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Old_Label : chr "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " ...
## $ Label : chr "BE B3 2/17 | G1 High1" "BE B3 2/17 | G1 High2" "BE B3 2/17 | G1 High3" "BE B3 2/17 | G1 High4" ...
## $ Population : Factor w/ 84 levels "High1","High13",..: 1 9 19 24 25 26 27 28 39 49 ...
## $ Genetic_Diversity : Factor w/ 3 levels "High","Medium",..: 1 1 1 1 1 1 1 3 3 3 ...
## $ Generation_Parents : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Generation_Eggs : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Date_Census : chr "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
## $ Date_Start_Eggs : chr "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
## $ Date_End_Eggs : chr "2/18/21" "2/18/21" "2/18/21" "2/18/21" ...
## $ Image_Box1 : chr "DSC_0155" "DSC_0156" "DSC_0157" "DSC_0158" ...
## $ Image_Box2 : chr NA NA NA NA ...
## $ MC_Box1 : num NA NA NA NA NA NA NA NA NA NA ...
## $ MC_Box2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ MC_Total_Adults : num 102 105 106 102 102 ...
## $ HC_Box1 : int NA NA NA NA NA NA NA NA NA NA ...
## $ HC_Box2 : int NA NA NA NA NA NA NA NA NA NA ...
## $ HC_Total_Adults : int 100 100 100 100 100 100 100 100 100 100 ...
## $ Nb_parents : num NA NA NA NA NA NA NA NA NA NA ...
## $ Growth_rate : num NA NA NA NA NA NA NA NA NA NA ...
## $ First_Throw : chr "no" "no" "no" "no" ...
## $ Extinction_finalstatus: Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ Less_than_5 : chr "no" "no" "no" "no" ...
## $ Nb_adults : int 100 100 100 100 100 100 100 100 100 100 ...
## $ Lambda : num NA NA NA NA NA NA NA NA NA NA ...
## $ obs : Factor w/ 504 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...head(data)## Block Old_Label Label Population Genetic_Diversity
## 1 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High1 High1 High
## 2 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High2 High2 High
## 3 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High3 High3 High
## 4 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High4 High4 High
## 5 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High5 High5 High
## 6 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High6 High6 High
## Generation_Parents Generation_Eggs Date_Census Date_Start_Eggs Date_End_Eggs
## 1 0 1 2/17/21 2/17/21 2/18/21
## 2 0 1 2/17/21 2/17/21 2/18/21
## 3 0 1 2/17/21 2/17/21 2/18/21
## 4 0 1 2/17/21 2/17/21 2/18/21
## 5 0 1 2/17/21 2/17/21 2/18/21
## 6 0 1 2/17/21 2/17/21 2/18/21
## Image_Box1 Image_Box2 MC_Box1 MC_Box2 MC_Total_Adults HC_Box1 HC_Box2
## 1 DSC_0155 <NA> NA NA 101.9 NA NA
## 2 DSC_0156 <NA> NA NA 104.7 NA NA
## 3 DSC_0157 <NA> NA NA 105.8 NA NA
## 4 DSC_0158 <NA> NA NA 101.9 NA NA
## 5 DSC_0159 <NA> NA NA 102.2 NA NA
## 6 DSC_0160 <NA> NA NA 101.7 NA NA
## HC_Total_Adults Nb_parents Growth_rate First_Throw Extinction_finalstatus
## 1 100 NA NA no no
## 2 100 NA NA no no
## 3 100 NA NA no no
## 4 100 NA NA no no
## 5 100 NA NA no no
## 6 100 NA NA no no
## Less_than_5 Nb_adults Lambda obs
## 1 no 100 NA 1
## 2 no 100 NA 2
## 3 no 100 NA 3
## 4 no 100 NA 4
## 5 no 100 NA 5
## 6 no 100 NA 6tail(data)## Block Old_Label Label Population
## 571 Block5 B5 - Backup Fam H 8/25 BE B5 | G6 Med 19 Med19
## 572 Block5 B5 - Backup Fam I 8/25 BE B5 | G6 Med 20 Med20
## 573 Block5 B5-Backup Fam A 8/25 BE B5 | G6 Med 21 Med21
## 574 Block5 B5-Backup Fam B 8/25 BE B5 | G6 Med 22 Med22
## 575 Block5 B5-Backup Fam C 8/25 BE B5 | G6 Med 23 Med23
## 576 Block5 B5-Backup Fam D 8/25 BE B5 | G6 Med 24 Med24
## Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 571 Medium 5 6 8/25/21
## 572 Medium 5 6 8/25/21
## 573 Medium 5 6 8/25/21
## 574 Medium 5 6 8/25/21
## 575 Medium 5 6 8/25/21
## 576 Medium 5 6 8/25/21
## Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2 MC_Box1 MC_Box2
## 571 8/25/21 8/26/21 DSC_0273 DSC_0274 2.00000 2.00000
## 572 8/25/21 8/26/21 <NA> <NA> NA NA
## 573 8/25/21 8/26/21 DSC_0277 DSC_0278 29.76978 32.92378
## 574 8/25/21 8/26/21 <NA> <NA> NA NA
## 575 8/25/21 8/26/21 DSC_0267 DSC_0268 25.31021 46.10092
## 576 8/25/21 8/26/21 DSC_0283 DSC_0284 71.09684 55.37443
## MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 571 4.0 2 2 4 8 0.500000
## 572 0.0 NA NA NA 0 NA
## 573 62.7 23 23 46 35 1.314286
## 574 0.0 NA NA NA NA NA
## 575 71.4 20 36 56 40 1.400000
## 576 126.5 76 57 133 41 3.243902
## First_Throw Extinction_finalstatus Less_than_5 Nb_adults Lambda obs
## 571 no no yes 4 0.500000 499
## 572 extinct yes yes NA NA 500
## 573 no no no 46 1.314286 501
## 574 extinct yes yes NA NA 502
## 575 no no no 56 1.400000 503
## 576 no no no 133 3.243902 5041.2 Remaining heterozygosity
#Upload growth rate end of experiment
data_he <- read.table(file=here::here("data", "He_ReMaining_BeatricePop.csv"), sep=",", header=TRUE)
data_he$Population <- as.factor(data_he$Population)
data_he <- data_he[,c(1,3)]
### Additional: He remaining
# New vector for the lost during exp
data_he$He_remain_exp <- NA
# He lost during exp
for(i in levels(data$Population)){
temp_data_pop <- subset(data, Population==i)
ifelse(sum(temp_data_pop$Generation_Parents)!=15,print("ERROR"),print(i))
temp_he_remaining_gen <- 1 - (1/(2*temp_data_pop$Nb_parents))
temp_he_remain_exp <- prod(temp_he_remaining_gen, na.rm = TRUE)
data_he$He_remain_exp[data_he$Population==i] <- temp_he_remain_exp
rm(temp_he_remain_exp,temp_he_remaining_gen,temp_data_pop)
}## [1] "High1"
## [1] "High13"
## [1] "High14"
## [1] "High15"
## [1] "High16"
## [1] "High17"
## [1] "High18"
## [1] "High19"
## [1] "High2"
## [1] "High20"
## [1] "High21"
## [1] "High22"
## [1] "High23"
## [1] "High24"
## [1] "High25"
## [1] "High26"
## [1] "High27"
## [1] "High28"
## [1] "High3"
## [1] "High32"
## [1] "High33"
## [1] "High34"
## [1] "High35"
## [1] "High4"
## [1] "High5"
## [1] "High6"
## [1] "High7"
## [1] "Low1"
## [1] "Low10"
## [1] "Low11"
## [1] "Low12"
## [1] "Low13"
## [1] "Low14"
## [1] "Low15"
## [1] "Low16"
## [1] "Low17"
## [1] "Low18"
## [1] "Low19"
## [1] "Low2"
## [1] "Low20"
## [1] "Low21"
## [1] "Low22"
## [1] "Low24"
## [1] "Low25"
## [1] "Low26"
## [1] "Low27"
## [1] "Low28"
## [1] "Low29"
## [1] "Low3"
## [1] "Low30"
## [1] "Low31"
## [1] "Low32"
## [1] "Low33"
## [1] "Low34"
## [1] "Low35"
## [1] "Low36"
## [1] "Low4"
## [1] "Low5"
## [1] "Low7"
## [1] "Low8"
## [1] "Low9"
## [1] "Med1"
## [1] "Med10"
## [1] "Med11"
## [1] "Med12"
## [1] "Med13"
## [1] "Med14"
## [1] "Med15"
## [1] "Med17"
## [1] "Med18"
## [1] "Med19"
## [1] "Med2"
## [1] "Med20"
## [1] "Med21"
## [1] "Med22"
## [1] "Med23"
## [1] "Med24"
## [1] "Med3"
## [1] "Med4"
## [1] "Med5"
## [1] "Med6"
## [1] "Med7"
## [1] "Med8"
## [1] "Med9"# Remaining He at the end of the experiment
data_he$He_remain_end <- data_he$He_remain * data_he$He_remain_exp1.3 Phenotyping dataset
#Upload growth rate end of experiment
data_phenotyping <- read.table(file=here::here("data", "Adaptation_Data.csv"), sep=",", header=TRUE)
#Factors variables
data_phenotyping$ID_Rep <- as.factor(data_phenotyping$ID_Rep)
data_phenotyping$Week <- as.factor(data_phenotyping$Week)
data_phenotyping$Block <- as.factor(data_phenotyping$Block)
data_phenotyping$Population <- as.factor(data_phenotyping$Population)
#Revalue
data_phenotyping$Genetic_Diversity <- as.factor(data_phenotyping$Divsersity)
data_phenotyping$Genetic_Diversity <- plyr::revalue(data_phenotyping$Genetic_Diversity, c("Med"="Medium"))
data_phenotyping$Genetic_Diversity <- factor(data_phenotyping$Genetic_Diversity, levels = c("High", "Medium", "Low"))
#Add sum total
data_phenotyping$Nb_ttx <- rowSums(data_phenotyping[,11:13], na.rm=TRUE)
#Proportion
data_phenotyping$p_adults <- data_phenotyping$Adults / data_phenotyping$Nb_ttx
data_phenotyping$p_pupae <- data_phenotyping$Pupae / data_phenotyping$Nb_ttx
data_phenotyping$p_larvae <- data_phenotyping$Larvae / data_phenotyping$Nb_ttx
data_phenotyping_all <- data_phenotyping
data_phenotyping <- data_phenotyping[data_phenotyping$Week=="5-week",]1.4 Merge dataset He
#Add heterozygosity
data <- merge(x = data, y = data_he, by = "Population", all.x=TRUE)
#Add heterozygosity
data_phenotyping <- merge(x = data_phenotyping, y = data_he, by = "Population", all.x=TRUE)1.5 Subset dataset generation
#Data beginning
data_G2 <- data[data$Generation_Eggs==2,]
#Data end
data_G6 <- data[data$Generation_Eggs==6,]1.6 Merge dataset phenotyping
## Merge dataset phenotyping
#merge dataset
temp_Start <- data_G2[,c("Block","Population",
"Genetic_Diversity","Lambda", "He_remain", "He_remain_end" )]
temp_Start$ID_Rep <- 1
temp_Start$Phenotyping <- "Initial"
temp_end <- data_phenotyping[,c("Block","Population",
"Genetic_Diversity","Lambda","ID_Rep", "He_remain", "He_remain_end" )]
temp_end$Phenotyping <- "Final"
data_evolution_Lambda <- rbind(temp_Start,temp_end)
data_evolution_Lambda$Phenotyping <- factor(data_evolution_Lambda$Phenotyping,
levels = c("Initial", "Final"))2 PLOT
2.1 Population size
PLOT_Pop_size<-ggplot2::ggplot(data, aes(x = factor(Generation_Eggs),
y = Nb_adults,
group = Population,
colour = Genetic_Diversity)) +
geom_point(position = position_dodge(0.4), size = 0.8, alpha = 0.7) +
geom_line(position = position_dodge(0.4), size = 0.2, alpha = 0.6) +
ylab("Population size") +
labs(color="Genetic diversity") +
xlab("Generation") +
theme_LO
PLOT_Pop_size
# cowplot::save_plot(file = here::here("figures","OldPlot_PopulationSize.pdf"),PLOT_Pop_size, base_height = 10/cm(1),
# base_width = 15/cm(1), dpi = 610)
SUM_Popsize <- Rmisc::summarySE(data,
measurevar = c("Nb_adults"),
groupvars = c("Generation_Eggs",
"Genetic_Diversity"),
na.rm=TRUE)
SUM_Popsize## Generation_Eggs Genetic_Diversity N Nb_adults sd se ci
## 1 1 High 27 100.00000 0.00000 0.000000 0.00000
## 2 1 Medium 23 100.00000 0.00000 0.000000 0.00000
## 3 1 Low 34 100.00000 0.00000 0.000000 0.00000
## 4 2 High 27 344.29630 73.77294 14.197610 29.18360
## 5 2 Medium 23 282.04348 100.75735 21.009360 43.57075
## 6 2 Low 34 282.61765 78.53006 13.467794 27.40043
## 7 3 High 27 151.85185 80.96899 15.582489 32.03027
## 8 3 Medium 23 118.73913 88.54491 18.462891 38.28969
## 9 3 Low 34 104.61765 85.13341 14.600260 29.70445
## 10 4 High 26 114.00000 80.20224 15.728954 32.39439
## 11 4 Medium 23 62.65217 64.50483 13.450188 27.89398
## 12 4 Low 34 46.38235 39.63241 6.796903 13.82840
## 13 5 High 27 103.62963 51.56486 9.923661 20.39838
## 14 5 Medium 19 64.68421 49.88437 11.444259 24.04350
## 15 5 Low 30 53.16667 45.98207 8.395139 17.16999
## 16 6 High 27 147.66667 42.60282 8.198916 16.85311
## 17 6 Medium 18 73.61111 44.05496 10.383855 21.90802
## 18 6 Low 27 79.88889 44.43520 8.551559 17.57798PLOT_Pop_size_average <- PLOT_Pop_size +
geom_point(data = SUM_Popsize, aes(x = factor(Generation_Eggs),
y = Nb_adults,
group = Genetic_Diversity,
colour = Genetic_Diversity),
size = 5, position = position_dodge(0.4)) +
geom_errorbar(data = SUM_Popsize, aes(x = factor(Generation_Eggs),
ymin = Nb_adults-ci, ymax = Nb_adults+ci,
group = Genetic_Diversity,
colour = Genetic_Diversity),
size = 1, width=.2, position = position_dodge(0.4)) +
geom_line(data = SUM_Popsize, aes(x = factor(Generation_Eggs),
y = Nb_adults,
group = Genetic_Diversity,
colour = Genetic_Diversity),
linetype = "dashed", size = 1.5, position = position_dodge(0.4))
PLOT_Pop_size_average
# cowplot::save_plot(file = here::here("figures","OldPlot_PopulationSize_average.pdf"), PLOT_Pop_size_average, base_height = 10/cm(1), base_width = 15/cm(1), dpi = 610)2.2 Extinction
## Extinction yes no
PLOT_Extinction<-ggplot2::ggplot(data, aes(x = factor(Generation_Eggs),
y = Nb_adults,
group = Population,
colour = Extinction_finalstatus)) +
facet_wrap(~Genetic_Diversity, ncol=3) +
#geom_point(position = position_dodge(0.4), size = 0.4, alpha = 0.7) +
geom_line(position = position_dodge(0.1), size = 0.4, alpha = 0.8) +
ylab("Population size") +
scale_color_manual(values = c("black", "red")) +
xlab("Generation") +
theme(legend.position = "none") +
theme_LO
PLOT_Extinction
# #
# cowplot::save_plot(here::here("figures","OldPlot_Extinction.pdf"), PLOT_Extinction, base_height = 6/cm(1),
# base_width = 20/cm(1), dpi = 610)
#
# #2.3 Growth rate
tapply(data_G2$Lambda,data_G2$Genetic_Diversity,mean)## High Medium Low
## 3.442963 2.820435 2.826176PLOT_Growth<-ggplot2::ggplot(data_G2, aes(x = Genetic_Diversity, y = Lambda,
colour = Genetic_Diversity)) +
geom_boxplot(aes(group = Genetic_Diversity),outlier.colour = NA) +
geom_jitter(width = 0.25, size = 1, alpha = 0.8) +
ylab("Growth rate") +
theme(legend.position = "none") +
xlab("Genetic diversity") +
theme_LO
PLOT_Growth
#
# cowplot::save_plot(here::here("figures","OldPlot_Growth_Start.pdf"), PLOT_Growth, base_height = 10/cm(1),
# base_width = 8/cm(1), dpi = 610)2.4 Growth rate evolution
SUM_MEAN_start_end <- Rmisc::summarySE(data_evolution_Lambda, measurevar = c("Lambda"),
groupvars = c("Genetic_Diversity","Phenotyping"),
na.rm=TRUE)
SUM_POP_start_end <- Rmisc::summarySE(data_evolution_Lambda, measurevar = c("Lambda"),
groupvars = c("Genetic_Diversity","Phenotyping",
"Population"),
na.rm=TRUE)
PLOT_Growth_startend<-ggplot(SUM_POP_start_end, aes(x = Phenotyping,
y = Lambda,
group = Population,
colour = Genetic_Diversity)) +
geom_point(position = position_dodge(0.4), size = 0.8, alpha = 0.7) +
geom_line(position = position_dodge(0.4), size = 0.2, alpha = 0.6) +
ylab("Growth rate") +
xlab("Phenotyping") +
theme_LO
PLOT_Growth_startend
PLOT_Growth_startend_average <- PLOT_Growth_startend +
geom_point(data = SUM_MEAN_start_end, aes(x = Phenotyping,
y = Lambda,
group = Genetic_Diversity,
colour = Genetic_Diversity),
size = 5, position = position_dodge(0.4)) +
geom_errorbar(data = SUM_MEAN_start_end, aes(x = Phenotyping,
ymin = Lambda-ci, ymax = Lambda+ci,
group = Genetic_Diversity,
colour = Genetic_Diversity),
size = 1, width=.2, position = position_dodge(0.4)) +
geom_line(data = SUM_MEAN_start_end, aes(x = Phenotyping,
y = Lambda,
group = Genetic_Diversity,
colour = Genetic_Diversity),
linetype = "dashed", size = 1.5, position = position_dodge(0.4))
PLOT_Growth_startend_average
#
# cowplot::save_plot(here::here("figures","OldPlot_Lambda_Evolution.pdf"),
# PLOT_Growth_startend_average, base_height = 10/cm(1),
# base_width = 15/cm(1), dpi = 610)2.5 Life stage proportion
SUM_Prop_adults<- Rmisc::summarySE(data_phenotyping_all,
measurevar = c("p_adults"),
groupvars = c("Genetic_Diversity","Week"),
na.rm = TRUE)
SUM_Prop_adults## Genetic_Diversity Week N p_adults sd se ci
## 1 High 4-week 50 0.6508098 0.17776618 0.025139934 0.05052059
## 2 High 5-week 105 0.9632040 0.05259208 0.005132461 0.01017786
## 3 Medium 4-week 24 0.7121270 0.15722597 0.032093616 0.06639070
## 4 Medium 5-week 49 0.9563136 0.06160975 0.008801393 0.01769639
## 5 Low 4-week 30 0.6568606 0.19029917 0.034743716 0.07105888
## 6 Low 5-week 59 0.9582021 0.07845154 0.010213520 0.02044458SUM_Prop_pupae<- Rmisc::summarySE(data_phenotyping_all,
measurevar = c("p_pupae"),
groupvars = c("Genetic_Diversity","Week"),
na.rm = TRUE)
SUM_Prop_pupae## Genetic_Diversity Week N p_pupae sd se ci
## 1 High 4-week 50 0.28825555 0.15629878 0.022103985 0.044419621
## 2 High 5-week 105 0.01458769 0.02407445 0.002349426 0.004658999
## 3 Medium 4-week 24 0.23019148 0.14506173 0.029610601 0.061254195
## 4 Medium 5-week 49 0.02283167 0.03485130 0.004978757 0.010010462
## 5 Low 4-week 30 0.29380227 0.17562943 0.032065401 0.065581108
## 6 Low 5-week 59 0.01675308 0.02459640 0.003202178 0.006409856SUM_Prop_larvae<- Rmisc::summarySE(data_phenotyping_all,
measurevar = c("p_larvae"),
groupvars = c("Genetic_Diversity","Week"),
na.rm = TRUE)
SUM_Prop_larvae## Genetic_Diversity Week N p_larvae sd se ci
## 1 High 4-week 50 0.06093466 0.04274850 0.006045551 0.012148989
## 2 High 5-week 105 0.02220828 0.04324043 0.004219833 0.008368088
## 3 Medium 4-week 24 0.05768150 0.05999020 0.012245449 0.025331640
## 4 Medium 5-week 49 0.02085474 0.04462923 0.006375605 0.012819012
## 5 Low 4-week 30 0.04933713 0.08198954 0.014969174 0.030615398
## 6 Low 5-week 59 0.02504479 0.06691319 0.008711355 0.017437671SUM_Prop <- data.frame(SUM_Prop_adults[,1:4],
p_pupae=SUM_Prop_pupae$p_pupae,
p_larvae=SUM_Prop_larvae$p_larvae)
SUM_Prop<-tidyr::gather(SUM_Prop,"Stage", "Proportion",4:6,-"Genetic_Diversity")
SUM_Prop## Genetic_Diversity Week N Stage Proportion
## 1 High 4-week 50 p_adults 0.65080979
## 2 High 5-week 105 p_adults 0.96320403
## 3 Medium 4-week 24 p_adults 0.71212702
## 4 Medium 5-week 49 p_adults 0.95631359
## 5 Low 4-week 30 p_adults 0.65686059
## 6 Low 5-week 59 p_adults 0.95820213
## 7 High 4-week 50 p_pupae 0.28825555
## 8 High 5-week 105 p_pupae 0.01458769
## 9 Medium 4-week 24 p_pupae 0.23019148
## 10 Medium 5-week 49 p_pupae 0.02283167
## 11 Low 4-week 30 p_pupae 0.29380227
## 12 Low 5-week 59 p_pupae 0.01675308
## 13 High 4-week 50 p_larvae 0.06093466
## 14 High 5-week 105 p_larvae 0.02220828
## 15 Medium 4-week 24 p_larvae 0.05768150
## 16 Medium 5-week 49 p_larvae 0.02085474
## 17 Low 4-week 30 p_larvae 0.04933713
## 18 Low 5-week 59 p_larvae 0.02504479SUM_Prop$Stage <- factor(SUM_Prop$Stage, levels = c("p_adults", "p_pupae", "p_larvae"))
PLOT_prop <- ggplot2::ggplot(data = SUM_Prop, aes(x = Genetic_Diversity,
y = Proportion, fill = Stage)) +
facet_wrap(~ Week, ncol=2) +
geom_bar(stat="identity") +
xlab("Level of genetic diversity") +
labs(color="Stage of individuals") +
ylab("Proportion of each stage") +
scale_fill_manual(values= c("#4EABB9", "#E95C2B", "#CBA938"),
breaks = c("p_adults", "p_pupae", "p_larvae"),
labels = c( "Adult", "Pupae","Larvae")) +
theme_LO
PLOT_prop
# cowplot::save_plot(here::here("figures","OldPlot_Life_stage_proportion.pdf"), PLOT_prop,
# base_height = 8/cm(1), base_width = 16/cm(1), dpi = 610)
# 2.6 Remaining heterozygosity
# All but with: He beginning for both and with pseudoreplication for final
ggplot2::ggplot(data_evolution_Lambda, aes(x = He_remain, y = Lambda,
colour = Genetic_Diversity)) +
facet_wrap(~ Phenotyping, ncol = 1)+
geom_point(size = 2, alpha = 0.8) +
ylab("Growth rate") +
theme(legend.position = "none") +
xlab("Remaining heterozygosity") +
theme_LO
PLOT_He_Initial <- ggplot2::ggplot(data_G2, aes(x = He_remain, y = Lambda,
colour = Genetic_Diversity)) +
geom_point(size = 2, alpha = 0.8) +
ylab("Growth rate") +
theme(legend.position = "none") +
xlab("Remaining heterozygosity") +
theme_LO
PLOT_He_Initial
# Remove pseudoreplication: average of each population
SUM_POP_He_Mean<- Rmisc::summarySE(data_evolution_Lambda[data_evolution_Lambda$Phenotyping=="Final",],
measurevar = c("Lambda"),
groupvars = c("Genetic_Diversity","He_remain_end",
"Population"),
na.rm=TRUE)
PLOT_He_Final <- ggplot2::ggplot(SUM_POP_He_Mean, aes(x = He_remain_end, y = Lambda,
colour = Genetic_Diversity)) +
geom_point(size = 2, alpha = 0.8) +
ylab("Growth rate") +
theme(legend.position = "none") +
xlab("Remaining heterozygosity") +
theme_LO
PLOT_He_Final
# ALL WITHOUT PSEUDO
SUM_POP_He_ALL<- Rmisc::summarySE(data_evolution_Lambda,
measurevar = c("Lambda"),
groupvars = c("Genetic_Diversity","Phenotyping", "He_remain", "He_remain_end",
"Population"),
na.rm=TRUE)
SUM_POP_He_ALL$HE <- ifelse(SUM_POP_He_ALL$Phenotyping=="Initial",SUM_POP_He_ALL$He_remain,SUM_POP_He_ALL$He_remain_end)
PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda,
colour = Genetic_Diversity,
shape = Phenotyping)) +
geom_point(size = 3, alpha = 0.8) +
scale_shape_manual(values = c(1, 19)) +
ylab("Growth rate") +
xlab("Remaining heterozygosity") +
theme_LO
PLOT_He_ALL
# PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",], aes(x = HE, y = Lambda,
# colour = Genetic_Diversity,
# shape = Phenotyping)) +
# geom_point(size = 3, alpha = 0.8) +
# scale_shape_manual(values = c(19)) +
# ylab("Growth rate") +
# xlab("Remaining heterozygosity") +
# theme_LO
# PLOT_He_ALL3 ANALYSIS
3.1 Probability of extinction
############
###### Clean dataset
############
#Prepare clean dataset
data_proba_extinction <- data_G6[,c("Block","Population","Genetic_Diversity","Extinction_finalstatus")]
data_proba_extinction$Extinction_finalstatus <- as.factor(data_proba_extinction$Extinction_finalstatus)
data_proba_extinction$y <- ifelse(data_proba_extinction$Extinction_finalstatus == "yes", 1, 0)
############
###### Analysis
############
#Analysis
mod0 <- lme4::glmer(y ~ Genetic_Diversity + Block + (1|Population),
data = data_proba_extinction, family = binomial)
summary(mod0)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: y ~ Genetic_Diversity + Block + (1 | Population)
## Data: data_proba_extinction
##
## AIC BIC logLik deviance df.resid
## 58.8 73.4 -23.4 46.8 78
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1010 -0.3208 -0.2215 0.0000 4.5137
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 0 0
## Number of obs: 84, groups: Population, 84
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -22.3274 5077.2816 -0.004 0.99649
## Genetic_DiversityMedium 19.3132 5077.2815 0.004 0.99696
## Genetic_DiversityLow 19.5192 5077.2815 0.004 0.99693
## BlockBlock4 0.7402 1.2714 0.582 0.56046
## BlockBlock5 3.0006 1.1265 2.664 0.00773 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Gnt_DM Gnt_DL BlckB4
## Gntc_DvrstM -1.000
## Gntc_DvrstL -1.000 1.000
## BlockBlock4 0.000 0.000 0.000
## BlockBlock5 0.000 0.000 0.000 0.737
## convergence code: 0
## boundary (singular) fit: see ?isSingular#Test genetic diversity effect
mod1 <- lme4::glmer(y ~ Block + (1|Population),
data = data_proba_extinction, family = binomial)
anova(mod0, mod1, test = "Chisq")## Data: data_proba_extinction
## Models:
## mod1: y ~ Block + (1 | Population)
## mod0: y ~ Genetic_Diversity + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 4 66.903 76.626 -29.451 58.903
## mod0 6 58.833 73.418 -23.417 46.833 12.069 2 0.002394 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#We keep genetic diversity
mod0 <- lme4::glmer(y ~ Genetic_Diversity-1 + Block + (1|Population),
data = data_proba_extinction, family = binomial)
summary(mod0)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: y ~ Genetic_Diversity - 1 + Block + (1 | Population)
## Data: data_proba_extinction
##
## AIC BIC logLik deviance df.resid
## 58.8 73.4 -23.4 46.8 78
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1010 -0.3208 -0.2215 0.0000 4.5137
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 5.55e-15 7.45e-08
## Number of obs: 84, groups: Population, 84
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## Genetic_DiversityHigh -3.737e+01 8.535e+06 0.000 1.00000
## Genetic_DiversityMedium -3.014e+00 1.129e+00 -2.669 0.00760 **
## Genetic_DiversityLow -2.808e+00 1.066e+00 -2.635 0.00843 **
## BlockBlock4 7.402e-01 1.271e+00 0.582 0.56046
## BlockBlock5 3.001e+00 1.127e+00 2.664 0.00773 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## Gnt_DH Gnt_DM Gnt_DL BlckB4
## Gntc_DvrstM 0.000
## Gntc_DvrstL 0.000 0.767
## BlockBlock4 0.000 -0.724 -0.786
## BlockBlock5 0.000 -0.843 -0.871 0.737
## convergence code: 0
## boundary (singular) fit: see ?isSingular############
###### Predicted value
############
#Extract probability of extinction
(p_high <- plogis(-3.950e+01))## [1] 7.004352e-18(p_med <- plogis(-3.014e+00))## [1] 0.04679739(p_low <- plogis(-2.808e+00))## [1] 0.0568934predict(mod0, type = "response",
re.form = NA,
newdata = data.frame(Genetic_Diversity = levels(data_proba_extinction$Genetic_Diversity),
Block = "Block4"))## 1 2 3
## 2.220446e-16 9.329490e-02 1.122446e-01############
###### Posthoc
############
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High -36.12 8534537 Inf -20378273 20378200
## Medium -1.77 1 Inf -3 0
## Low -1.56 1 Inf -3 0
##
## Results are averaged over the levels of: Block
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium -34.352 8534537 Inf 0.000 1.0000
## High - Low -34.558 8534537 Inf 0.000 1.0000
## Medium - Low -0.206 1 Inf -0.274 0.9594
##
## Results are averaged over the levels of: Block
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates3.2 Population size G6
############
###### Analysis
############
#Analysis
# mod0 <- lme4::glmer(Nb_adults ~ Genetic_Diversity + Block + (1|Population),
# data = data_G6, family = "poisson")
mod0 <- lm(log(Nb_adults) ~ Genetic_Diversity + Block,
data = data_G6[data_G6$Extinction_finalstatus == "no",])
mod1 <- lm(log(Nb_adults) ~ Block ,
data = data_G6[data_G6$Extinction_finalstatus == "no",])
anova(mod0, mod1) # Effect of genetic diversity ## Analysis of Variance Table
##
## Model 1: log(Nb_adults) ~ Genetic_Diversity + Block
## Model 2: log(Nb_adults) ~ Block
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 66 29.706
## 2 68 42.118 -2 -12.412 13.789 9.914e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod0 <- lm(log(Nb_adults) ~ Genetic_Diversity + Block,
data = data_G6[data_G6$Extinction_finalstatus == "no",])
summary(mod0)##
## Call:
## lm(formula = log(Nb_adults) ~ Genetic_Diversity + Block, data = data_G6[data_G6$Extinction_finalstatus ==
## "no", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.2417 -0.2947 0.1088 0.4381 1.2624
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.11532 0.17980 28.451 < 2e-16 ***
## Genetic_DiversityMedium -0.94813 0.20540 -4.616 1.86e-05 ***
## Genetic_DiversityLow -0.80532 0.18719 -4.302 5.72e-05 ***
## BlockBlock4 -0.02077 0.18447 -0.113 0.9107
## BlockBlock5 -0.53922 0.21414 -2.518 0.0142 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6709 on 66 degrees of freedom
## Multiple R-squared: 0.3362, Adjusted R-squared: 0.2959
## F-statistic: 8.356 on 4 and 66 DF, p-value: 1.623e-05############
###### Posthoc
############
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df lower.CL upper.CL
## High 4.93 0.130 66 4.61 5.25
## Medium 3.98 0.159 66 3.59 4.37
## Low 4.12 0.136 66 3.79 4.46
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df t.ratio p.value
## High - Medium 0.948 0.205 66 4.616 0.0001
## High - Low 0.805 0.187 66 4.302 0.0002
## Medium - Low -0.143 0.207 66 -0.690 0.7704
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates3.3 Population size over time
PLOT_Pop_size_average
fit <- glm(Nb_adults ~ Generation_Eggs, data = data, family = "poisson")
#second degree
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE), data = data, family = "poisson")
#third degree
fit3 <- glm(Nb_adults~stats::poly(Generation_Eggs,3,raw=TRUE), data = data, family = "poisson")
#fourth degree
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE), data = data, family = "poisson")
#fourth degree
fit5 <- glm(Nb_adults~stats::poly(Generation_Eggs,5,raw=TRUE), data = data, family = "poisson")
#exp
fit6 <- glm(Nb_adults~exp(Generation_Eggs), data = data, family = "poisson")
#log
fit7 <- glm(Nb_adults~log(Generation_Eggs), data = data, family = "poisson")
#generate range of 50 numbers starting from 30 and ending at 160
xx <- seq(1,6, length=100)
plot(data$Generation_Eggs,
data$Nb_adults,pch=19)
lines(xx, predict(fit, data.frame(Generation_Eggs=xx),type = "response"), col="red")
lines(xx, predict(fit2, data.frame(Generation_Eggs=xx),type = "response"), col="green")
lines(xx, predict(fit3, data.frame(Generation_Eggs=xx),type = "response"), col="blue")
lines(xx, predict(fit4, data.frame(Generation_Eggs=xx),type = "response"), col="purple")
lines(xx, predict(fit5, data.frame(Generation_Eggs=xx),type = "response"), col="yellow")
lines(xx, predict(fit6, data.frame(Generation_Eggs=xx),type = "response"), col="orange")
lines(xx, predict(fit7, data.frame(Generation_Eggs=xx),type = "response"), col="pink")
# Check the r squared
summary(fit)$adj.r.squared## NULLsummary(fit2)$adj.r.squared## NULLsummary(fit3)$adj.r.squared## NULLsummary(fit4)$adj.r.squared # best r square## NULLsummary(fit5)$adj.r.squared ## NULLsummary(fit6)$adj.r.squared ## NULLsummary(fit7)$adj.r.squared ## NULL#
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE), data = data, family = "poisson")
#If we dont consider Gen=1
PLOT_Pop_size_average
fit <- lm(log(Nb_adults+1) ~ Generation_Eggs, data = data[data$Generation_Eggs>=2,])
#second degree
fit2 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,2,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#third degree
fit3 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,3,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#fourth degree
fit4 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,4,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#fourth degree
fit5 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,5,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#exp
fit6 <- lm(log(Nb_adults+1)~exp(Generation_Eggs), data = data[data$Generation_Eggs>=2,])
#log
fit7 <- lm(log(Nb_adults+1)~log(Generation_Eggs), data = data[data$Generation_Eggs>=2,])
#generate range of 50 numbers starting from 30 and ending at 160
xx <- seq(1,6, length=100)
plot(data[data$Generation_Eggs>=2,]$Generation_Eggs,
log(data[data$Generation_Eggs>=2,]$Nb_adults+1),pch=19,ylim=c(0,8))
lines(xx, predict(fit, data.frame(Generation_Eggs=xx)), col="red")
lines(xx, predict(fit2, data.frame(Generation_Eggs=xx)), col="green")
lines(xx, predict(fit3, data.frame(Generation_Eggs=xx)), col="blue")
lines(xx, predict(fit4, data.frame(Generation_Eggs=xx)), col="purple")
lines(xx, predict(fit5, data.frame(Generation_Eggs=xx)), col="yellow")
lines(xx, predict(fit6, data.frame(Generation_Eggs=xx)), col="orange")
lines(xx, predict(fit7, data.frame(Generation_Eggs=xx)), col="pink")
# Check the r squared
summary(fit)$adj.r.squared## [1] 0.1281606summary(fit2)$adj.r.squared## [1] 0.3033469summary(fit3)$adj.r.squared## [1] 0.3017434summary(fit4)$adj.r.squared # best r square## [1] 0.3031656summary(fit5)$adj.r.squared ## [1] 0.3031656summary(fit6)$adj.r.squared ## [1] 0.01432146summary(fit7)$adj.r.squared ## [1] 0.1763784# data$gen_square <- data$Generation_Eggs^2
# fit2 <- lm(log(Nb_adults+1)~poly(Generation_Eggs,2,raw=TRUE), data = data[data$Generation_Eggs>=2,])
# fit2 <- lm(log(Nb_adults+1)~ Generation_Eggs + gen_square, data = data[data$Generation_Eggs>=2,])
# #Add the interaction with Genetic diversity
# fit2 <- glm(Nb_adults~ Generation_Eggs*Genetic_Diversity +
# gen_square*Genetic_Diversity,
# data = data, family = "poisson")
#
#
# summary(fit2)
#######################
####################### MODEL INCLUDING GENETIC DIVERSITY
#######################
#Test oversdispersion
fit4 <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + (1|Block), data = data, family = "poisson")
blmeco::dispersion_glmer(fit4)## [1] 5.69289mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity +
(1|Block) + (1|obs), data = data, family = "poisson")
blmeco::dispersion_glmer(mod)## [1] 1.075746#Convergence issue
max(abs(with(mod@optinfo$derivs, solve(Hessian, gradient)))) #should be <0.001## [1] 0.07273844#Predict data
filldata <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity),100),
Generation_Eggs = rep(seq(1,6, length = 100),each=3),
Block = "Block4",
Estimates = NA)
filldata$Estimates <- predict(fit4, newdata=filldata, type = "response")
#Plot
PLOT_Pop_size <- ggplot2::ggplot(data = data) +
geom_line(data = filldata, aes(x = Generation_Eggs,
y = Estimates,
colour = Genetic_Diversity), size = 1.3) +
geom_point(data = data,
aes(x = Generation_Eggs,
y = Nb_adults,
colour = Genetic_Diversity),
position = position_dodge(0.4), size = 1.4, alpha = 0.6) +
ylab("Population size") +
labs(color="Genetic diversity") +
xlab("Generation") +
theme_LO
PLOT_Pop_size
# fit2 <- lm(log(Nb_adults+1)~ Generation_Eggs*Genetic_Diversity +
# gen_square*Genetic_Diversity,
# data = data[data$Generation_Eggs>=2,])
# fit2_test <- lm(log(Nb_adults+1)~ Generation_Eggs+Genetic_Diversity +
# gen_square,
# data = data[data$Generation_Eggs>=2,])
mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity +
(1|Block) + (1|obs), data = data, family = "poisson")
mod1 <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)+Genetic_Diversity +
(1|Block) + (1|obs), data = data, family = "poisson")
anova(mod,mod1, test = "Chisq")## Data: data
## Models:
## mod1: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity +
## mod1: (1 | Block) + (1 | obs)
## mod: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +
## mod: (1 | Block) + (1 | obs)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 9 5541.3 5579 -2761.7 5523.3
## mod 17 5527.0 5598 -2746.5 5493.0 30.365 8 0.0001822 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Signifcant interaction
summary(mod)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +
## (1 | Block) + (1 | obs)
## Data: data
##
## AIC BIC logLik deviance df.resid
## 5527.0 5598.0 -2746.5 5493.0 466
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.12632 -0.04460 0.00513 0.06208 0.29142
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 0.58420 0.7643
## Block (Intercept) 0.07837 0.2799
## Number of obs: 483, groups: obs, 483; Block, 3
##
## Fixed effects:
## Estimate
## (Intercept) -1.987774
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 10.850818
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -5.143118
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 0.941456
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -0.058946
## Genetic_DiversityMedium -0.598092
## Genetic_DiversityLow -0.917922
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 1.230996
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium -0.776495
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.157397
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium -0.010610
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow 1.729870
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow -0.976354
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow 0.171344
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow -0.009432
## Std. Error
## (Intercept) 1.224862
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 1.974563
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 1.029881
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 0.213080
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 0.015165
## Genetic_DiversityMedium 1.822946
## Genetic_DiversityLow 1.616291
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 2.977179
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 1.559898
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.324231
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.023172
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow 2.630789
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow 1.374871
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow 0.285225
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow 0.020358
## z value
## (Intercept) -1.623
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 5.495
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -4.994
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 4.418
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -3.887
## Genetic_DiversityMedium -0.328
## Genetic_DiversityLow -0.568
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 0.413
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium -0.498
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.485
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium -0.458
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow 0.658
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow -0.710
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow 0.601
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow -0.463
## Pr(>|z|)
## (Intercept) 0.104620
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 3.90e-08
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 5.92e-07
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 9.95e-06
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 0.000101
## Genetic_DiversityMedium 0.742843
## Genetic_DiversityLow 0.570090
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 0.679257
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 0.618635
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.627359
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.647046
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow 0.510829
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow 0.477616
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow 0.548017
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow 0.643123
##
## (Intercept)
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 ***
## Genetic_DiversityMedium
## Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 1.3357 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?summary(mod1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity +
## (1 | Block) + (1 | obs)
## Data: data
##
## AIC BIC logLik deviance df.resid
## 5541.3 5579.0 -2761.7 5523.3 474
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.09399 -0.05085 0.01382 0.07441 0.23786
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 0.62340 0.7896
## Block (Intercept) 0.07624 0.2761
## Number of obs: 483, groups: obs, 483; Block, 3
##
## Fixed effects:
## Estimate Std. Error z value
## (Intercept) -2.164709 0.781262 -2.771
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 11.963215 1.252559 9.551
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -5.794296 0.658426 -8.800
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 1.063434 0.137051 7.759
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -0.066353 0.009797 -6.773
## Genetic_DiversityMedium -0.510577 0.095339 -5.355
## Genetic_DiversityLow -0.641010 0.086094 -7.445
## Pr(>|z|)
## (Intercept) 0.00559 **
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 < 2e-16 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 < 2e-16 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 8.53e-15 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 1.26e-11 ***
## Genetic_DiversityMedium 8.54e-08 ***
## Genetic_DiversityLow 9.66e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) s::(G_E,4,r=TRUE)1 s::(G_E,4,r=TRUE)2
## s::(G_E,4,r=TRUE)1 -0.964
## s::(G_E,4,r=TRUE)2 0.940 -0.993
## s::(G_E,4,r=TRUE)3 -0.914 0.976 -0.995
## s::(G_E,4,r=TRUE)4 0.888 -0.957 0.984
## Gntc_DvrstM -0.063 0.009 -0.010
## Gntc_DvrstL -0.066 0.005 -0.006
## s::(G_E,4,r=TRUE)3 s::(G_E,4,r=TRUE)4 Gnt_DM
## s::(G_E,4,r=TRUE)1
## s::(G_E,4,r=TRUE)2
## s::(G_E,4,r=TRUE)3
## s::(G_E,4,r=TRUE)4 -0.997
## Gntc_DvrstM 0.011 -0.011
## Gntc_DvrstL 0.007 -0.007 0.491
## convergence code: 0
## Model failed to converge with max|grad| = 0.338298 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?emmeans::emmeans(mod, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High 4.57 0.203 Inf 4.08 5.05
## Medium 3.94 0.212 Inf 3.44 4.45
## Low 3.71 0.197 Inf 3.24 4.18
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium 0.624 0.182 Inf 3.434 0.0017
## High - Low 0.854 0.163 Inf 5.229 <.0001
## Medium - Low 0.231 0.177 Inf 1.306 0.3916
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates#######################
####################### MODEL INCLUDING ONLY RESCUED POPULATIONS
#######################
#Same but only for population not extinct
fit_fin <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity,
data = data[data$Extinction_finalstatus=="no",], family = "poisson")
filldata$Estimates_Withoutextpop <- predict(fit_fin, newdata=filldata, type = "response")
PLOT_Pop_size <- ggplot2::ggplot(data = data) +
geom_line(data = filldata, aes(x = Generation_Eggs,
y = Estimates_Withoutextpop,
colour = Genetic_Diversity), size = 1.5, linetype = "longdash") +
geom_point(data = data[data$Extinction_finalstatus=="no",],
aes(x = Generation_Eggs,
y = Nb_adults,
colour = Genetic_Diversity),
position = position_dodge(0.4), size = 1.4, alpha = 0.5) +
ylab("Population size") +
labs(color="Genetic diversity") +
xlab("Generation") +
theme_LO
PLOT_Pop_size
mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity +
(1|Block) + (1|obs),
data = data[data$Extinction_finalstatus=="no",], family = "poisson")
mod1 <- lme4::glmer(Nb_adults~poly(Generation_Eggs,4,raw=TRUE)+Genetic_Diversity +
(1|Block) + (1|obs),
data = data[data$Extinction_finalstatus=="no",], family = "poisson")
anova(mod,mod1, test = "Chisq")## Data: data[data$Extinction_finalstatus == "no", ]
## Models:
## mod1: Nb_adults ~ poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity +
## mod1: (1 | Block) + (1 | obs)
## mod: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +
## mod: (1 | Block) + (1 | obs)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 9 4759.9 4796.3 -2370.9 4741.9
## mod 17 4756.3 4825.1 -2361.1 4722.3 19.599 8 0.01196 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mod)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +
## (1 | Block) + (1 | obs)
## Data: data[data$Extinction_finalstatus == "no", ]
##
## AIC BIC logLik deviance df.resid
## 4756.3 4825.1 -2361.1 4722.3 407
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.57063 -0.04879 0.00393 0.06967 0.29650
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 0.36624 0.6052
## Block (Intercept) 0.02456 0.1567
## Number of obs: 424, groups: obs, 424; Block, 3
##
## Fixed effects:
## Estimate
## (Intercept) -1.9250178
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 10.7535583
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -5.0922588
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 0.9316535
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -0.0583188
## Genetic_DiversityMedium 1.1892488
## Genetic_DiversityLow 0.2794991
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium -2.0341925
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 1.0112012
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium -0.2039827
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.0137179
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow -0.4034287
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow 0.0777534
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow -0.0101264
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow 0.0006044
## Std. Error
## (Intercept) 0.9705855
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 1.5722049
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 0.8204596
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 0.1698215
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 0.0120893
## Genetic_DiversityMedium 1.5049889
## Genetic_DiversityLow 1.3816159
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 2.4448787
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 1.2751857
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.2640143
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.0188078
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow 2.2495506
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow 1.1751948
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow 0.2435505
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow 0.0173582
## z value
## (Intercept) -1.983
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 6.840
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -6.207
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 5.486
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -4.824
## Genetic_DiversityMedium 0.790
## Genetic_DiversityLow 0.202
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium -0.832
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 0.793
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium -0.773
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.729
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow -0.179
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow 0.066
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow -0.042
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow 0.035
## Pr(>|z|)
## (Intercept) 0.0473
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 7.93e-12
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 5.41e-10
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 4.11e-08
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 1.41e-06
## Genetic_DiversityMedium 0.4294
## Genetic_DiversityLow 0.8397
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 0.4054
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 0.4278
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.4397
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.4658
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow 0.8577
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow 0.9472
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow 0.9668
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow 0.9722
##
## (Intercept) *
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 ***
## Genetic_DiversityMedium
## Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 3.25364 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?emmeans::emmeans(mod, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High 4.53 0.133 Inf 4.21 4.85
## Medium 4.30 0.151 Inf 3.94 4.66
## Low 4.01 0.137 Inf 3.68 4.33
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium 0.230 0.153 Inf 1.507 0.2874
## High - Low 0.523 0.140 Inf 3.746 0.0005
## Medium - Low 0.293 0.157 Inf 1.861 0.1503
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates#Signifcant interaction 3.4 Evolution of fitness over time
############
###### Distribution
############
hist(data_evolution_Lambda$Lambda, breaks=50)
# Square root
msqrt <- lme4::lmer(sqrt(Lambda) ~ Genetic_Diversity*Phenotyping +
Block + (1|Population),
data = data_evolution_Lambda)
summary(msqrt)## Linear mixed model fit by REML ['lmerMod']
## Formula: sqrt(Lambda) ~ Genetic_Diversity * Phenotyping + Block + (1 |
## Population)
## Data: data_evolution_Lambda
##
## REML criterion at convergence: 252.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4979 -0.4737 0.0672 0.5620 2.4161
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 0.0295 0.1718
## Residual 0.1040 0.3225
## Number of obs: 301, groups: Population, 88
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.779696 0.081683 21.788
## Genetic_DiversityMedium -0.171161 0.103522 -1.653
## Genetic_DiversityLow -0.150003 0.094035 -1.595
## PhenotypingFinal -0.243939 0.070714 -3.450
## BlockBlock4 0.102764 0.063883 1.609
## BlockBlock5 0.003821 0.070206 0.054
## Genetic_DiversityMedium:PhenotypingFinal -0.251050 0.110515 -2.272
## Genetic_DiversityLow:PhenotypingFinal -0.154999 0.101256 -1.531
##
## Correlation of Fixed Effects:
## (Intr) Gnt_DM Gnt_DL PhntyF BlckB4 BlckB5 G_DM:P
## Gntc_DvrstM -0.600
## Gntc_DvrstL -0.655 0.506
## PhntypngFnl -0.685 0.534 0.588
## BlockBlock4 -0.462 0.054 0.041 0.028
## BlockBlock5 -0.424 0.006 0.009 0.002 0.488
## Gntc_DvM:PF 0.428 -0.745 -0.377 -0.641 -0.028 0.050
## Gntc_DvL:PF 0.444 -0.373 -0.733 -0.699 0.012 0.088 0.453# log
mlog <- lme4::lmer(log(Lambda+1) ~ Genetic_Diversity*Phenotyping +
Block + (1|Population),
data=data_evolution_Lambda)
# Normal
mnorm <- lme4::lmer(Lambda ~ Genetic_Diversity*Phenotyping +
Block + (1|Population),
data=data_evolution_Lambda)
## Compare AIC
AIC(mlog, msqrt,mnorm)## df AIC
## mlog 10 207.0029
## msqrt 10 272.0671
## mnorm 10 848.9685#Square root a better fits to the data
#But we compare different values: Growth rate and Growth rate+1
#Distribution Residuals
hist(residuals(mlog),col="yellow",freq=F)
hist(residuals(msqrt),col="yellow",freq=F) #Seems the better fit
hist(residuals(mnorm),col="yellow",freq=F)
#Test normality
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mlog), "norm"))
g1$kstest ## 1-mle-norm
## "rejected"g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(msqrt), "norm"))
g1$kstest ## 1-mle-norm
## "not rejected"g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mnorm), "norm"))
g1$kstest## 1-mle-norm
## "not rejected"# QQ plot
# plot(mlog, which = 2)
# plot(msqrt, which = 2) #Seems the better fit!
# plot(mnorm, which = 2) #Seems the better fit but very high AIC
qqnorm(resid(mlog))
qqline(resid(mlog))
qqnorm(resid(msqrt))
qqline(resid(msqrt))
qqnorm(resid(mnorm))
qqline(resid(mnorm)) #best fit
############
###### Model selection
############
m0 <- lme4::lmer(Lambda ~ Genetic_Diversity*Phenotyping +
Block + (1|Population) ,
data=data_evolution_Lambda)
m1 <- lme4::lmer(Lambda ~ Genetic_Diversity+Phenotyping +
Block + (1|Population) ,
data=data_evolution_Lambda)
anova(m0,m1,test="Chisq") #Remove interaction## Data: data_evolution_Lambda
## Models:
## m1: Lambda ~ Genetic_Diversity + Phenotyping + Block + (1 | Population)
## m0: Lambda ~ Genetic_Diversity * Phenotyping + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 8 832.93 862.59 -408.46 816.93
## m0 10 834.25 871.32 -407.12 814.25 2.6798 2 0.2619#Same result with:
# m1 <- lmer(Lambda ~ Genetic_Diversity+Phenotyping +
# (1|Block) + (1|Population) + (1|ID_Rep:Population),
# data=data_evolution_Lambda)
# OR
# m1 <- lmer(Lambda ~ Genetic_Diversity+Phenotyping +
# (1|Block) + (1|Population) ,
# data=data_evolution_Lambda)
m2 <- lme4::lmer(Lambda ~ Phenotyping +
Block + (1|Population),
data=data_evolution_Lambda)
m3 <- lme4::lmer(Lambda ~ Genetic_Diversity +
Block + (1|Population),
data=data_evolution_Lambda)
anova(m1,m2,test="Chisq") #Keep phenotyping## Data: data_evolution_Lambda
## Models:
## m2: Lambda ~ Phenotyping + Block + (1 | Population)
## m1: Lambda ~ Genetic_Diversity + Phenotyping + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m2 6 855.30 877.54 -421.65 843.30
## m1 8 832.93 862.59 -408.46 816.93 26.369 2 1.88e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m1,m3,test="Chisq") #Keep genetic diversity## Data: data_evolution_Lambda
## Models:
## m3: Lambda ~ Genetic_Diversity + Block + (1 | Population)
## m1: Lambda ~ Genetic_Diversity + Phenotyping + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m3 7 896.49 922.44 -441.24 882.49
## m1 8 832.93 862.59 -408.46 816.93 65.558 1 5.642e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1############
###### Posthoc
############
mod_final <- lme4::lmer(Lambda ~ Genetic_Diversity+Phenotyping +
Block + (1|Population),
data=data_evolution_Lambda)
emmeans::emmeans(mod_final, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df lower.CL upper.CL
## High 3.05 0.120 69.4 2.76 3.34
## Medium 2.14 0.143 87.6 1.79 2.48
## Low 2.32 0.123 99.8 2.02 2.62
##
## Results are averaged over the levels of: Phenotyping, Block
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df t.ratio p.value
## High - Medium 0.916 0.184 72.8 4.988 <.0001
## High - Low 0.734 0.170 78.2 4.316 0.0001
## Medium - Low -0.182 0.187 87.5 -0.972 0.5963
##
## Results are averaged over the levels of: Phenotyping, Block
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmeans::emmeans(mod_final, list(pairwise ~ Phenotyping), adjust = "tukey") ### $`emmeans of Phenotyping`
## Phenotyping emmean SE df lower.CL upper.CL
## Initial 3.00 0.1072 254 2.76 3.24
## Final 2.01 0.0837 109 1.82 2.20
##
## Results are averaged over the levels of: Genetic_Diversity, Block
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 2 estimates
##
## $`pairwise differences of Phenotyping`
## contrast estimate SE df t.ratio p.value
## Initial - Final 0.993 0.117 250 8.474 <.0001
##
## Results are averaged over the levels of: Genetic_Diversity, Block
## Degrees-of-freedom method: kenward-roger#3.5 Fitness G6
m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Population) + (1|ID_Rep:Population),
family = "poisson", data = data_phenotyping)
summary(m0)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
## Data: data_phenotyping
##
## AIC BIC logLik deviance df.resid
## 2154.4 2171.3 -1072.2 2144.4 212
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1441 -0.1230 0.0254 0.1360 0.3923
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_Rep:Population (Intercept) 0.2041 0.4518
## Population (Intercept) 0.2807 0.5298
## Number of obs: 217, groups: ID_Rep:Population, 217; Population, 77
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.2792 0.1105 38.715 < 2e-16 ***
## Genetic_DiversityMedium -0.7386 0.1812 -4.076 4.59e-05 ***
## Genetic_DiversityLow -0.5250 0.1664 -3.156 0.0016 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Gnt_DM
## Gntc_DvrstM -0.611
## Gntc_DvrstL -0.667 0.413# Equivalent to:
# m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Population/ID_Rep),
# family = "poisson", data = data_5week)
#dispersion_lme4::glmer(m0) #dispersion ok
#Note: (1|School) + (1|Class:School) is equivalent to (1|School/Class)
m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Population) + (1|ID_Rep:Population),
family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1 + (1|Population) + (1|ID_Rep:Population),
family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") # difference## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Population) + (1 | ID_Rep:Population)
## m0: Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 3 2167.5 2177.6 -1080.7 2161.5
## m0 5 2154.4 2171.3 -1072.2 2144.4 17.04 2 0.0001994 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m0)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
## Data: data_phenotyping
##
## AIC BIC logLik deviance df.resid
## 2154.4 2171.3 -1072.2 2144.4 212
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1441 -0.1230 0.0254 0.1360 0.3923
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_Rep:Population (Intercept) 0.2041 0.4518
## Population (Intercept) 0.2807 0.5298
## Number of obs: 217, groups: ID_Rep:Population, 217; Population, 77
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.2792 0.1105 38.715 < 2e-16 ***
## Genetic_DiversityMedium -0.7386 0.1812 -4.076 4.59e-05 ***
## Genetic_DiversityLow -0.5250 0.1664 -3.156 0.0016 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Gnt_DM
## Gntc_DvrstM -0.611
## Gntc_DvrstL -0.667 0.413############
###### Posthoc
############
emmeans::emmeans(m0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High 4.28 0.111 Inf 4.02 4.54
## Medium 3.54 0.143 Inf 3.20 3.88
## Low 3.75 0.124 Inf 3.46 4.05
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium 0.739 0.181 Inf 4.076 0.0001
## High - Low 0.525 0.166 Inf 3.156 0.0046
## Medium - Low -0.214 0.189 Inf -1.132 0.4943
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates#test double nested
m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Genetic_Diversity/Population/ID_Rep),
family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1 + (1|Genetic_Diversity/Population/ID_Rep),
family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") ## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Genetic_Diversity/Population/ID_Rep)
## m0: Adults ~ Genetic_Diversity + (1 | Genetic_Diversity/Population/ID_Rep)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 4 2160.9 2174.4 -1076.4 2152.9
## m0 6 2156.4 2176.7 -1072.2 2144.4 8.4441 2 0.01467 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Genetic_Diversity/ID_Rep),
family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1 + (1|Genetic_Diversity/ID_Rep),
family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") ## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Genetic_Diversity/ID_Rep)
## m0: Adults ~ Genetic_Diversity + (1 | Genetic_Diversity/ID_Rep)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 3 2189.9 2200.0 -1091.9 2183.9
## m0 5 2183.0 2199.9 -1086.5 2173.0 10.858 2 0.004387 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.6 Fitness G2
m0 <- glm(Nb_adults ~ Genetic_Diversity + Block, family = "poisson", data = data_G2)
summary(m0)##
## Call:
## glm(formula = Nb_adults ~ Genetic_Diversity + Block, family = "poisson",
## data = data_G2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -18.7057 -2.7621 0.2126 3.3983 11.4454
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.73890 0.01497 383.312 < 2e-16 ***
## Genetic_DiversityMedium -0.19408 0.01628 -11.920 < 2e-16 ***
## Genetic_DiversityLow -0.19278 0.01460 -13.205 < 2e-16 ***
## BlockBlock4 0.10767 0.01579 6.817 9.3e-12 ***
## BlockBlock5 0.17743 0.01600 11.088 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2379.6 on 83 degrees of freedom
## Residual deviance: 2027.9 on 79 degrees of freedom
## AIC: 2667.2
##
## Number of Fisher Scoring iterations: 4m1<- glm(Nb_adults ~ Block, family = "poisson", data = data_G2)
anova(m0,m1,test="Chisq") # difference## Analysis of Deviance Table
##
## Model 1: Nb_adults ~ Genetic_Diversity + Block
## Model 2: Nb_adults ~ Block
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 79 2027.9
## 2 81 2241.4 -2 -213.52 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1############
###### Posthoc
############
emmeans::emmeans(m0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High 5.834 0.01051 Inf 5.809 5.859
## Medium 5.640 0.01243 Inf 5.610 5.670
## Low 5.641 0.01022 Inf 5.617 5.666
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium 0.1941 0.0163 Inf 11.920 <.0001
## High - Low 0.1928 0.0146 Inf 13.205 <.0001
## Medium - Low -0.0013 0.0161 Inf -0.081 0.9964
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates3.7 Remaining He difference
Rmisc::summarySE(data_G2,
measurevar = c("He_remain"),
groupvars = c("Genetic_Diversity"),
na.rm=TRUE)## Genetic_Diversity N He_remain sd se ci
## 1 High 27 0.9986008 0.000000000 0.000000000 0.000000000
## 2 Medium 23 0.6791246 0.006059307 0.001263453 0.002620241
## 3 Low 34 0.5441887 0.012687849 0.002175948 0.004427000mfull <- lm(He_remain ~ Genetic_Diversity, data=data_G2)
mless <- lm(He_remain ~ 1, data=data_G2)
anova(mfull,mless) #Keep effect## Analysis of Variance Table
##
## Model 1: He_remain ~ Genetic_Diversity
## Model 2: He_remain ~ 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 81 0.0061
## 2 83 3.1868 -2 -3.1807 21048 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1############
###### Posthoc
############
emmeans::emmeans(mfull, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df lower.CL upper.CL
## High 0.9986 0.001673 81 0.9945 1.0027
## Medium 0.6791 0.001812 81 0.6747 0.6835
## Low 0.5442 0.001491 81 0.5406 0.5478
##
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df t.ratio p.value
## High - Medium 0.319 0.00247 81 129.527 <.0001
## High - Low 0.454 0.00224 81 202.800 <.0001
## Medium - Low 0.135 0.00235 81 57.498 <.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimatesRmisc::summarySE(data_G2,
measurevar = c("He_remain"),
groupvars = c("Genetic_Diversity"),
na.rm=TRUE)## Genetic_Diversity N He_remain sd se ci
## 1 High 27 0.9986008 0.000000000 0.000000000 0.000000000
## 2 Medium 23 0.6791246 0.006059307 0.001263453 0.002620241
## 3 Low 34 0.5441887 0.012687849 0.002175948 0.004427000#Calcul for end:
Rmisc::summarySE(data_G2[data_G2$Extinction_finalstatus!="yes",],
measurevar = c("He_remain_end"),
groupvars = c("Genetic_Diversity"),
na.rm=TRUE)## Genetic_Diversity N He_remain_end sd se ci
## 1 High 27 0.9697741 0.01868272 0.003595491 0.007390637
## 2 Medium 18 0.6482803 0.02453741 0.005783522 0.012202164
## 3 Low 26 0.5146307 0.02759467 0.005411760 0.011145728mfull <- lm(He_remain_end ~ Genetic_Diversity, data=data_G2[data_G2$Extinction_finalstatus!="yes",])
mless <- lm(He_remain_end ~ 1, data=data_G2[data_G2$Extinction_finalstatus!="yes",])
anova(mfull,mless) #Keep effect## Analysis of Variance Table
##
## Model 1: He_remain_end ~ Genetic_Diversity
## Model 2: He_remain_end ~ 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 68 0.03835
## 2 70 2.91180 -2 -2.8735 2547.7 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1############
###### Posthoc
############
emmeans::emmeans(mfull, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df lower.CL upper.CL
## High 0.970 0.00457 68 0.959 0.981
## Medium 0.648 0.00560 68 0.635 0.662
## Low 0.515 0.00466 68 0.503 0.526
##
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df t.ratio p.value
## High - Medium 0.321 0.00723 68 44.491 <.0001
## High - Low 0.455 0.00653 68 69.754 <.0001
## Medium - Low 0.134 0.00728 68 18.355 <.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimates3.8 Growth rate and He
#
# #fit first degree polynomial equation:
# fit <- lm(Lambda ~ HE, data = SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",])
# #same fit as: fit5 <- lm(Lambda~HE+HE^2, data = SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",])
# #same fit as: fit6 <- lm(Lambda~HE^2, data = SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",])
#
# #second degree
# fit2 <- lm(Lambda~poly(HE,2,raw=TRUE), data = SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",])
# #third degree
# fit3 <- lm(Lambda~poly(HE,3,raw=TRUE), data = SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",])
# #fourth degree
# fit4 <- lm(Lambda~poly(HE,4,raw=TRUE), data = SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",])
# #Other fit
# fit7 <- lm(Lambda~exp(HE), data = SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",])
# fit8 <- lm(Lambda~log(HE), data = SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",])
#
# #generate range of 50 numbers starting from 30 and ending at 160
# xx <- seq(0.3,1, length=100)
# plot(SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",]$HE,
# SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",]$Lambda,pch=19,ylim=c(0,5))
# lines(xx, predict(fit, data.frame(HE=xx)), col="red")
# lines(xx, predict(fit2, data.frame(HE=xx)), col="green")
# lines(xx, predict(fit3, data.frame(HE=xx)), col="blue")
# lines(xx, predict(fit4, data.frame(HE=xx)), col="purple")
# # lines(xx, predict(fit5, data.frame(HE=xx)), col="black")
# # lines(xx, predict(fit6, data.frame(HE=xx)), col="yellow")
# lines(xx, predict(fit7, data.frame(HE=xx)), col="brown")
# lines(xx, predict(fit8, data.frame(HE=xx)), col="orange")
#